---
title: "Moral Values Across Countries"
output:
flexdashboard::flex_dashboard:
orientation: columns
vertical_layout: fill
source_code: embed
theme: united
---
```{r setup, include=FALSE}
# global options
knitr::opts_chunk$set(echo = FALSE,
tidy = TRUE,
cache = FALSE,
message = FALSE,
error = FALSE,
warning = FALSE)
# packages
library(flexdashboard)
library(here)
library(rio)
library(tidyverse)
library(magrittr)
library(lme4)
library(lmerTest)
library(colorBlindness)
library(dotwhisker)
library(tidytext)
library(ggeffects)
library(see)
theme_set(theme_minimal()) # set theme
options(scipen=999) # remove scientific notation
```
```{r wrangling, include = FALSE}
# import data
df <- import(here("data", "ALL_MFQ30.csv"), # moral values, countries, & sex
setclass = "tbl_df") %>%
janitor::clean_names()
df_c <- import(here("data", "Data_S1_sec.csv"), # country-level variables
setclass = "tbl_df") %>%
janitor::clean_names()
# data wrangling
df %<>%
drop_na() %>%
mutate(
across(where(is.double), as.numeric),
across(where(is.character), as.factor),
sex = recode(sex,
`1` = "Male",
`0` = "Female",
.default = NA_character_),
indiv = rowMeans(
select(., harm_avg, fairness_avg) # individualizing moral foundations
),
bind = rowMeans(
select(., ingroup_avg:purity_avg) # binding moral foundations
)
)
# check data structure and variables
str(df)
# descriptive statistics by country
c_desc <-
df %>%
pivot_longer(cols = c(indiv, bind),
names_to = "vars",
values_to = "val"
) %>%
select(country, vars, val) %>%
group_by(country, vars) %>%
summarise(mean = mean(val, na.rm = TRUE),
sd = sd(val, na.rm = TRUE),
min = min(val, na.rm = TRUE),
max = max(val, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(vars = fct_recode(vars,
Individualizing = "indiv",
Binding = "bind"
)
)
# descriptive statistics by country and sex
c_s_desc <-
df %>%
filter(country != "Poland") %>% # Poland has missing data in sex.
pivot_longer(cols = c(indiv, bind),
names_to = "vars",
values_to = "val"
) %>%
group_by(country, sex, vars) %>%
summarise(mean = mean(val, na.rm = TRUE),
sd = sd(val, na.rm = TRUE),
min = min(val, na.rm = TRUE),
max = max(val, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(vars = fct_recode(vars,
Individualizing = "indiv",
Binding = "bind"
)
)
```
# Values X Country
Sidebar {.sidebar}
-----------------------------------------------------------------------
**Data**
Data used in this dashboard come from the second study of [Atari et al. (2020)](http://dx.doi.org/10.1098/rspb.2020.1201). It's a publicly available dataset, which can be downloaded from [Kaggle](https://www.kaggle.com/tunguz/sex-differences-in-moral-judgements-67-countries).
This study has data on moral values in 19 countries. There is also a second dataset with countr-level variables (see below for the list of the variables).
**Summary of the Visualizations**
The first set of plots (on this page) represents the average scores for individualizing and binding foundations across countries.
The second and third set of plots visualize the effect of sex and several country-level variables on the aforementioned foundations.
*Country-level variables*
Population Sex Ratio
Individualism
Masculinity
Gender Equality
Human Development Index
Overall Life Satisfaction Index
Column {data-width=600}
-----------------------------------------------------------------------
### Final Version
```{r}
gmeans <-
c_desc %>%
group_by(vars) %>%
summarise(m = mean(mean))
c_desc %>%
ggplot() +
geom_vline(data = gmeans,
aes(xintercept = m),
linetype = 2,
alpha = .6) +
geom_col(
aes(mean, reorder_within(country, mean, vars),
fill = country,
alpha = .9
)
) +
scale_y_reordered() +
scale_x_continuous(expand = c(0, 0)
) +
facet_wrap(~vars,
scales = "free_y",
ncol = 2) +
theme(
plot.title.position = "plot",
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
legend.position = "none",
axis.text.y = element_text(color = "black",
size = 11),
axis.text.x = element_text(color = "black",
size = 9),
axis.title = element_blank()
) +
labs(
title = "Endorsement of Individualizing and Binding Moral Values Across Countries",
caption = "Vertical lines represent the average of all countries."
)
#Colorizing each country may be is fine, but it kind of defeats its own purpose when there are over 20 individualized colors... I believe grouping countries based on their regions (e.g. Europe, Asia, North America, and etc) will improve interpretability of your plot by a lot! In Lab 3, we used colors to highlight the region of USA each state is from.
```
Column {data-width=400}
-----------------------------------------------------------------------
### Initial version
```{r}
c_desc %>%
ggplot() +
geom_col(
aes(mean, country, fill = vars),
position = "dodge"
) +
theme(
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()
) +
scale_x_continuous(expand = c(0, 0))
```
> This is the first version of the plot. There is a lot to work on here. Both axes and legend labels seem confusing. X-axis scale is also not complete. Also, it is hard to see any patterns without using facet_wrap and sorting the values.
### Revised version
```{r}
c_desc %>%
ggplot() +
geom_col(
aes(mean, reorder_within(country, mean, vars)
)
) +
scale_y_reordered() +
scale_x_continuous(expand = c(0, 0)
) +
facet_wrap(~vars,
scales = "free_y",
ncol = 2) +
theme(
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.y = element_text(color = "black",
size = 11),
axis.text.x = element_text(color = "black",
size = 9),
axis.title = element_blank()
) +
labs(
title = "Endorsement of Individualizing and Binding Moral Values Across Countries"
)
```
> The revised version looks much better. You can see interesting patterns such as Spain's scores. The colors are terrible, though. Also, it'd be nice to see the grand mean to have a general reference category.
# Predictors of moral values
Sidebar {.sidebar}
-----------------------------------------------------------------------
**Model**
I conducted a multilevel regression for each of the moral values where sex was the Level 1 predictor, the country-level variables were the Level 2 predictors, country had a random intercept, and a random slope for sex was used.
*R code*:
lmer(indiv ~ sex +
pop_sex_ratio + individualism +
masculinity + gender_eqality +
human_development_index +
overall_life_satisfaction_index +
(sex | country),
data = df)
**Interpretation**
The plot describes the fixed effects for the models described above. The coefficients are unstandardized and represented by the dots. The lines represent the 95% confidence intervals.
Looking at the plot, we can see that sex, gender equality, and overall life satisfaction are the significant predictors of individualizing moral values; whereas, only gender equality is a significant predictor of binding moral values.
Column {data-width=600}
-----------------------------------------------------------------------
```{r MLM, include=FALSE}
# merge country-level data with the main dataset
df <- left_join(df,
select(df_c, country:overall_life_satisfaction_index),
by = "country")
# MLM
model_i <- lmer(indiv ~ sex + # level 1 predictor
# level 2 predictors:
pop_sex_ratio + individualism + masculinity +
gender_eqality + human_development_index +
overall_life_satisfaction_index +
(sex | country), # random slope for sex
# random intercept for country
data = df
)
model_b <- lmer(bind ~ sex + # level 1 predictor
# level 2 predictors:
pop_sex_ratio + individualism + masculinity +
gender_eqality + human_development_index +
overall_life_satisfaction_index +
(sex | country), # random slope for sex
# random intercept for country
data = df
)
sjPlot::tab_model(model_i)
sjPlot::tab_model(model_b)
# extract coefficients
m_i_fixed <-
broom.mixed::tidy(model_i) %>%
filter(effect == "fixed",
term != "(Intercept)") %>%
select(-c(effect, group)) %>%
mutate(
term = recode(term,
`sexMale` = "Sex",
`pop_sex_ratio` = "Population Sex Ratio",
`individualism` = "Individualism",
`masculinity` = "Masculinity",
`gender_eqality` = "Gender Equality",
`human_development_index` = "Human Development Index",
`overall_life_satisfaction_index` = "Overall Life Satisfaction Index"),
model = "Individualizing",
) %>%
relocate(model, term)
m_b_fixed <-
broom.mixed::tidy(model_b) %>%
filter(effect == "fixed",
term != "(Intercept)") %>%
select(-c(effect, group)) %>%
mutate(
term = recode(term,
`sexMale` = "Sex",
`pop_sex_ratio` = "Population Sex Ratio",
`individualism` = "Individualism",
`masculinity` = "Masculinity",
`gender_eqality` = "Gender Equality",
`human_development_index` = "Human Development Index",
`overall_life_satisfaction_index` = "Overall Life Satisfaction Index"),
model = "Binding",
) %>%
relocate(model, term)
both_ms <- bind_rows(m_i_fixed, m_b_fixed)
```
### Final Version
```{r}
dwplot(both_ms,
dot_args = list(size = 2)
) +
ggtitle("Predicting moral values by sex and country-level predictors") +
xlab("Unstandardized Coefficient") +
geom_vline(xintercept = 0,
colour = "grey60",
linetype = 2) +
theme(plot.title = element_text(face = "bold", vjust = 3),
plot.title.position = "plot",
axis.text.y = element_text(color = "black",
size = 11),
legend.justification = c(0, 0),
legend.position = c(.65, .85),
legend.background = element_rect(colour = "grey80"),
legend.title = element_blank(),
panel.grid.major.y = element_blank()
) +
scale_x_continuous(n.breaks = 10) +
scale_color_manual(values = c("cornflowerblue", "#F8766D"))
#Great job! That said, I had an error with dw_plot() command for some reason... dwplot() works fine so I made some adjustments accordingly.
#I assume this visualizations cater to scholaristic audiences, so it would also be helpful to add APA style summaries of the regression analyses in another tab.
```
Column {data-width=400}
-----------------------------------------------------------------------
### Initial version
```{r}
dwplot(both_ms) +
ggtitle("Predicting moral values") +
xlab("Unstandardized Coefficient") +
geom_vline(xintercept = 0,
colour = "grey60",
linetype = 2)
```
> The initial dot-whisker plot for the fixed effects. Added a vertical line, which made the plot a bit easier to interpret. However, legend looks awful. Modifying the x axis should also help.
### Revised version
```{r}
dwplot(both_ms,
dot_args = list(size = 2)
) +
ggtitle("Predicting moral values by sex and country-level predictors") +
xlab("Unstandardized Coefficient") +
geom_vline(xintercept = 0,
colour = "grey60",
linetype = 2) +
theme(plot.title = element_text(face = "bold"),
plot.title.position = "plot",
axis.text.y = element_text(color = "black",
size = 11),
legend.justification = c(0, 0),
legend.position = c(.74, .85),
legend.background = element_rect(colour = "grey80"),
legend.title.align = .5
) +
scale_x_continuous(n.breaks = 10) +
scale_color_manual(values = c("blue", "red"))
```
> Easier to see the points. The legend makes sense now, but still could be better. Colors can be improved. Also, there are too many grid lines.
# Predictors of moral values for each country
Sidebar {.sidebar}
-----------------------------------------------------------------------
On this page, I visualize the predicted values of individualizing and binding moral values. Based on the earlier analysis conducted in the "Predictors of moral values" section, I only focus gender equality as the predictor. In fact, I reran the model as follows:
lmer(indiv ~ gender_eqality +
(sex | country),
data = df)
Column {data-width=600}
-----------------------------------------------------------------------
### Final Version
```{r}
# run the models
model1 <- lmer(indiv ~ gender_eqality + (sex|country), data = df)
model2 <- lmer(bind ~ gender_eqality + (sex|country), data = df)
# extract the predicted values
predicted1 <-
ggpredict(model1,
terms = c("gender_eqality", "country"),
type = "re")
predicted2 <-
ggpredict(model2,
terms = c("gender_eqality", "country"),
type = "re")
p11 <-
predicted1 %>%
ggplot(aes(x, predicted, color = group)) +
geom_line(size = 1) +
labs(x = "Gender Equality",
y = "Individualizing",
color = "Country") +
theme(axis.text = element_text(size = 10,
colour = "black")
) +
gghighlight::gghighlight(group %in% c("Poland", "Netherlands", "Hungary"))
p21 <-
predicted2 %>%
ggplot(aes(x, predicted, color = group)) +
geom_line(size = 1) +
labs(x = "Gender Equality",
y = "Binding",
color = "Country") +
theme(axis.text = element_text(size = 10,
colour = "black")
) +
gghighlight::gghighlight(group %in% c("Poland", "Netherlands", "Spain"))
ggpubr::ggarrange(p11, p21,
common.legend = TRUE,
legend = "bottom") %>%
ggpubr::annotate_figure(
top = ggpubr::text_grob("Countries at the top, middle, and bottom")
)
#The slope for the right side seems steeper, but I think it could be because the minimum and maximum values of y axis are set different for these tables. What about using the same values (min 2 and max 4, 0.5 increment each) for both tables?
#Similarly to the suggestion I made for prior plots, but could you group_by() and summarize() averages of these countries into regions (e.g. Europe, Asia, and etc)? That way we need to worry less about which countries to hightlight!
```
Column {data-width=400}
-----------------------------------------------------------------------
### Initial version
```{r}
# plot
p1 <-
predicted1 %>%
ggplot(aes(x, predicted, color = group)) +
geom_line(size = 1)
p2 <-
predicted2 %>%
ggplot(aes(x, predicted, color = group)) +
geom_line(size = 1)
ggpubr::ggarrange(p1, p2,
common.legend = TRUE,
legend = "bottom"
)
```
> The initial attempt to visualize the predicted values of moral foundations where the predictor is gender equality. I chose this predictor because it was the only that was significant for both outcomes. I am glad that it worked, but it needs improvement.
### Revised version
```{r}
p1 <-
p1 +
labs(x = "Gender Equality",
y = "Individualizing",
color = "Country") +
theme(axis.text = element_text(size = 10,
colour = "black"))
p2 <-
p2 +
labs(x = "Gender Equality",
y = "Binding",
color = "Country") +
theme(axis.text = element_text(size = 10,
colour = "black"))
ggpubr::ggarrange(p1, p2,
common.legend = TRUE,
legend = "bottom"
)
```
> This version is better, but it would be nicer to see which countries are at the top, middle, and bottom of the plot.
